home *** CD-ROM | disk | FTP | other *** search
/ Aminet 5 / Aminet 5 - March 1995.iso / Aminet / util / libs / type1beta5.lha / type1 / src / arith.c < prev    next >
C/C++ Source or Header  |  1994-12-17  |  7KB  |  249 lines
  1. /* $XConsortium: arith.c,v 1.4 94/03/22 19:08:54 gildea Exp $ */
  2. /* Copyright International Business Machines, Corp. 1991
  3.  * All Rights Reserved
  4.  * Copyright Lexmark International, Inc. 1991
  5.  * All Rights Reserved
  6.  *
  7.  * License to use, copy, modify, and distribute this software and its
  8.  * documentation for any purpose and without fee is hereby granted,
  9.  * provided that the above copyright notice appear in all copies and that
  10.  * both that copyright notice and this permission notice appear in
  11.  * supporting documentation, and that the name of IBM or Lexmark not be
  12.  * used in advertising or publicity pertaining to distribution of the
  13.  * software without specific, written prior permission.
  14.  *
  15.  * IBM AND LEXMARK PROVIDE THIS SOFTWARE "AS IS", WITHOUT ANY WARRANTIES OF
  16.  * ANY KIND, EITHER EXPRESS OR IMPLIED, INCLUDING, BUT NOT LIMITED TO ANY
  17.  * IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE,
  18.  * AND NONINFRINGEMENT OF THIRD PARTY RIGHTS.  THE ENTIRE RISK AS TO THE
  19.  * QUALITY AND PERFORMANCE OF THE SOFTWARE, INCLUDING ANY DUTY TO SUPPORT
  20.  * OR MAINTAIN, BELONGS TO THE LICENSEE.  SHOULD ANY PORTION OF THE
  21.  * SOFTWARE PROVE DEFECTIVE, THE LICENSEE (NOT IBM OR LEXMARK) ASSUMES THE
  22.  * ENTIRE COST OF ALL SERVICING, REPAIR AND CORRECTION.  IN NO EVENT SHALL
  23.  * IBM OR LEXMARK BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
  24.  * DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
  25.  * PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
  26.  * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
  27.  * THIS SOFTWARE.
  28.  */
  29.  
  30. /*
  31.  * ARITH Module - Portable Module for Multiple Precision Fixed Point Arithmetic
  32.  *
  33.  * This module provides division and multiplication of 64-bit fixed point
  34.  * numbers.  (To be more precise, the module works on numbers that take
  35.  * two 'longs' to store.  That is almost always equivalent to saying 64-bit
  36.  * numbers.)
  37.  *
  38.  * Note: it is frequently easy and desirable to recode these functions in
  39.  * assembly language for the particular processor being used, because
  40.  * assembly language, unlike C, will have 64-bit multiply products and
  41.  * 64-bit dividends.  This module is offered as a portable version.
  42.  *
  43.  * &author. Jeffrey B. Lotspiech (lotspiech@almaden.ibm.com) and Sten F. Andler
  44.  *
  45.  *
  46.  *
  47.  * Reference for all algorithms:  Donald E. Knuth, "The Art of Computer
  48.  * Programming, Volume 2, Semi-Numerical Algorithms," Addison-Wesley Co.,
  49.  * Massachusetts, 1969, pp. 229-279.
  50.  *
  51.  * Knuth talks about a 'digit' being an arbitrary sized unit and a number
  52.  * being a sequence of digits.  We'll take a digit to be a 'short'.
  53.  *
  54.  * The following assumption must be valid for these algorithms to work:
  55.  *    A 'long' is two 'short's.
  56.  *
  57.  * The following code is INDEPENDENT of:
  58.  *    The actual size of a short.
  59.  *    Whether shorts and longs are stored most significant byte
  60.  *          first or least significant byte first.
  61.  */
  62.  
  63.  
  64. /*
  65.  * Include Files
  66.  *
  67.  * The included files are:
  68.  */
  69. #ifndef T1GST
  70. #include "global.h"
  71. #endif
  72. #include <stdio.h>
  73.  
  74.  
  75. /**
  76.  * Structure Definitions
  77.  **/
  78. typedef struct
  79. {
  80.     long high;
  81.     unsigned long low;
  82. }
  83. doublelong;
  84.  
  85.  
  86. /*
  87.  * SHORTSIZE is the number of bits in a short; LONGSIZE is the number of
  88.  * bits in a long; MAXSHORT is the maximum unsigned short:
  89.  */
  90. #define     SHORTSIZE         (sizeof(short)*8)
  91. #define     LONGSIZE          (SHORTSIZE*2)
  92. #define     MAXSHORT          ((1<<SHORTSIZE)-1)
  93.  
  94.  
  95. /*
  96.  * ASSEMBLE concatenates two shorts to form a long:
  97.  */
  98. #define     ASSEMBLE(hi,lo)   ((((unsigned long)hi)<<SHORTSIZE)+(lo))
  99.  
  100.  
  101. /*
  102.  * HIGHDIGIT extracts the most significant short from a long; LOWDIGIT
  103.  * extracts the least significant short from a long:
  104.  */
  105. #define     HIGHDIGIT(u)      ((u)>>SHORTSIZE)
  106. #define     LOWDIGIT(u)       ((u)&MAXSHORT)
  107.  
  108.  
  109. /*
  110.  * SIGNBITON tests the high order bit of a long 'w':
  111.  */
  112. #define    SIGNBITON(w)   (((long)w)<0)
  113.  
  114.  
  115. /**
  116.  * Local prototypes
  117.  **/
  118. static void DLmult(doublelong *product, unsigned long u, unsigned long v);
  119.  
  120. /*
  121.  * DLrightshift() - Macro to Shift Double Long Right by N
  122.  */
  123. #define  DLrightshift(dl,N)  { \
  124.        dl.low = (dl.low >> N) + (((unsigned long) dl.high) << (LONGSIZE - N)); \
  125.        dl.high >>= N; \
  126. }
  127.  
  128.  
  129. /*
  130.  * Double Long Arithmetic
  131.  *
  132.  * DLmult() - Multiply Two Longs to Yield a Double Long
  133.  *
  134.  * The two multiplicands must be positive.
  135.  */
  136. //
  137. //OPTIMIZATION BY AMISH 4/26/93
  138. //In my test cases, this function was always called with
  139. //u OR v == 1 (NOT TOFRACTPEL(1), but rather a 32 bit #
  140. //with the first bit on. - (important distinction)
  141. //
  142. //So, calling DLmult is silly in these cases, since we
  143. //can just set product.low = u or v (whichever is not 1)
  144. //and product.high = 0.
  145. //
  146. //NOTE: This resulted in a barely measurable speedup
  147. //
  148. static void DLmult(doublelong *product, unsigned long u, unsigned long v)
  149. {
  150.     unsigned long u1, u2;    /* the digits of u */
  151.     unsigned long v1, v2;    /* the digits of v */
  152.     unsigned int w1, w2, w3, w4;    /* the digits of w */
  153.     unsigned long t;    /* temporary variable */
  154.  
  155. //THE FOLLOWING TWO IF STATEMENTS ADDED BY AMISH 4/26/93
  156.     if (u == 1)
  157.     {
  158.         product->high = 0;
  159.         product->low = v;
  160.         return;
  161.     }
  162.  
  163.     if (v == 1)
  164.     {
  165.         product->high = 0;
  166.         product->low = u;
  167.         return;
  168.     }
  169.  
  170.     u1 = HIGHDIGIT(u);
  171.     u2 = LOWDIGIT(u);
  172.     v1 = HIGHDIGIT(v);
  173.     v2 = LOWDIGIT(v);
  174.  
  175.     if (v2 == 0)
  176.         w4 = w3 = w2 = 0;
  177.     else
  178.     {
  179.         t = u2 * v2;
  180.         w4 = LOWDIGIT(t);
  181.         t = u1 * v2 + HIGHDIGIT(t);
  182.         w3 = LOWDIGIT(t);
  183.         w2 = HIGHDIGIT(t);
  184.     }
  185.  
  186.     if (v1 == 0)
  187.         w1 = 0;
  188.     else
  189.     {
  190.         t = u2 * v1 + w3;
  191.         w3 = LOWDIGIT(t);
  192.         t = u1 * v1 + w2 + HIGHDIGIT(t);
  193.         w2 = LOWDIGIT(t);
  194.         w1 = HIGHDIGIT(t);
  195.     }
  196.  
  197.     product->high = ASSEMBLE(w1, w2);
  198.     product->low = ASSEMBLE(w3, w4);
  199. }
  200.  
  201.  
  202. /*
  203.  * Fractional Pel Arithmetic
  204.  */
  205.  
  206. /*
  207.  * FPmult() - Multiply Two Fractional Pel Values
  208.  *
  209.  * This funtion first calculates w = u * v to "doublelong" precision.
  210.  * It then shifts w right by FRACTBITS bits, and checks that no
  211.  * overflow will occur when the resulting value is passed back as
  212.  * a fractpel.
  213.  */
  214. fractpel FPmult(fractpel u, fractpel v)
  215. {
  216.     doublelong w;
  217.     int negative = FALSE;    /* sign flag */
  218.  
  219.     if ((u == 0) || (v == 0))
  220.         return (0);
  221.  
  222.     if (u < 0)
  223.     {
  224.         u = -u;
  225.         negative = TRUE;
  226.     }
  227.     if (v < 0)
  228.     {
  229.         v = -v;
  230.         negative = !negative;
  231.     }
  232.  
  233.     if (u == TOFRACTPEL(1))
  234.         return ((negative) ? -v : v);
  235.     if (v == TOFRACTPEL(1))
  236.         return ((negative) ? -u : u);
  237.  
  238.     DLmult(&w, u, v);
  239.     DLrightshift(w, FRACTBITS);
  240.     if (w.high != 0 || SIGNBITON(w.low))
  241.     {
  242. //        IfTrace2(TRUE, "FPmult: overflow, %px%p\n", u, v);
  243.         w.low = TOFRACTPEL(MAXSHORT);
  244.     }
  245.  
  246.     return (fractpel)((negative) ? -w.low : w.low);
  247. }
  248.  
  249.